@doc """
A collection of useful mathematical functions by Julia language. Author: https://gitee.com/iamzhtr

In REPL use `?MathKits.f` get help on the function `f`, for example, `?MathKits.pnpoly!`. 

Available functions: `crossproduct`, `distance_to_segment`, `norm`/`norm!`, `norm2`/`norm2!`, `pnpoly`/`pnpoly!`, `product`/`product!`, `root_on_segment`, `between`, `betweeneq`
""" ->
module MathKits
using Markdown

# ----------------------------------- #
function between(x::T, a::T, b::T) where T <: Real
    return x > a && x < b
end

function betweeneq(x::T, a::T, b::T) where T <: Real
    return x >= a && x <= b
end

@doc """
`ans = crossproduct(a::Vector{T} where T <: Number, b::Vector{T} where T <: Number)`

Return the cross product of two vectors.
""" ->
function crossproduct(a::Vector{T} where T <: Number, b::Vector{T} where T <: Number)
    return [a[2]*b[3]-b[2]*a[3],a[3]*b[1]-a[1]*b[3],a[1]*b[2]-b[1]*a[2]]
end

@doc """
`ans = norm(v::Array{T} where T <: Number)`

Return the L2 norm of the vector consisting of all elements of the array `v`.
""" ->
norm(v::Array{T} where T <: Number) = sqrt(norm2!(v))

@doc """
`ans = norm!(v::Array{T} where T <: Number)`

Return the L2 norm of the vector consisting of all elements of the array `v`.
""" ->
norm!(v::Array{T} where T <: Number) = sqrt(norm2!(v))

@doc """
`ans = norm2(v::Array{T} where T <: Number)`

Return the sum of all squares of element of the array `v`, i.e., `norm2!(v) = norm!(v)^2`.
""" ->
norm2(v::Array{T} where T <: Number) = sum(v .* v)

@doc """
`ans = norm2!(v::Array{T} where T <: Number)`

Return the sum of all squares of element of the array `v`, i.e., `norm2!(v) = norm!(v)^2`.
""" ->
norm2!(v::Array{T} where T <: Number) = sum(v .* v)

@doc """
```julia
ans = pnpoly(vertices_x::Vector, vertices_y::Vector, point_x::Real, point_y::Real)
```
Find if a point lies within a 2D polygon. Return 0 if it lies inside, else return 1. This algorithm was proposed by W. Randolph Franklin: https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html

`vertices_x`/`vertices_y`: Vector of `x`/`y` of the polygon vertices. Note that `x[end] != x[1]` (very important).
 
`point_x`/`point_y`: `x`/`y` of the point.
""" ->  
function pnpoly(vertices_x::Vector, vertices_y::Vector, point_x::Real, point_y::Real)
    @assert length(vertices_x) == length(vertices_y)
    nvert = length(vertices_x)
    VertX = [vertices_x; vertices_x[1]]
    VertY = [vertices_y; vertices_y[1]]
    if ( (point_x < minimum(VertX)) || (point_x > maximum(VertX)) || (point_y < minimum(VertY)) || (point_y > maximum(VertY)) )
        c = 1
    else
        c = 1
        for i = 1:nvert
            j = i + 1
            if ( (VertY[i] > point_y) ⊻ (VertY[j] > point_y ) ) && ( point_x < (VertX[j] - VertX[i]) * (point_y - VertY[i]) / (VertY[j] - VertY[i]) + VertX[i] )
                # ⊻ : \xor
                c = 1 - c
            end
        end
    end
    # c = 0: inside, 1: outside
    return c
end

@doc """
```julia
ans = pnpoly!(vertices_x::Vector, vertices_y::Vector, point_x::Real, point_y::Real)
```
Find if a point lies within a 2D polygon. Return 0 if it lies inside, else return 1. This algorithm was proposed by W. Randolph Franklin: https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html

`vertices_x`/`vertices_y`: Vector of `x`/`y` of the polygon vertices. Note that `x[end] != x[1]` (very important).
 
`point_x`/`point_y`: `x`/`y` of the point.
""" ->  
function pnpoly!(vertices_x::Vector, vertices_y::Vector, point_x::Real, point_y::Real)
    @assert length(vertices_x) == length(vertices_y)
    nvert = length(vertices_x)
    VertX = [vertices_x; vertices_x[1]]
    VertY = [vertices_y; vertices_y[1]]
    if ( (point_x < minimum(VertX)) || (point_x > maximum(VertX)) || (point_y < minimum(VertY)) || (point_y > maximum(VertY)) )
        c = 1
    else
        c = 1
        for i = 1:nvert
            j = i + 1
            if ( (VertY[i] > point_y) ⊻ (VertY[j] > point_y ) ) && ( point_x < (VertX[j] - VertX[i]) * (point_y - VertY[i]) / (VertY[j] - VertY[i]) + VertX[i] )
                # ⊻ : \xor
                c = 1 - c
            end
        end
    end
    # c = 0: inside, 1: outside
    return c
end

@doc """
`ans = product(v::Array{T} where T <: Number)`

Return the product of all elements of the array `v`.
""" ->
function product(v::Array{T} where T <: Number)
    s = 1.
    for a in v
        s *= a
    end
    s
end

@doc """
`ans = product!(v::Array{T} where T <: Number)`

Return the product of all elements of the array `v`.
""" ->
function product!(v::Array{T} where T <: Number)
    s = 1.
    for a in v
        s *= a
    end
    s
end

@doc """
`root, lambda = root_on_segment(vertex1::Array{T} where T <: Number, vertex2::Array{T} where T <: Number, point::Array{T} where T <: Number)`

Return the `root` point and length ratio `lambda` of `point` on the segment from `vertex1` to `vertex2`. The length ratio 'lambda = ||root - vertex1|| / ||vertex2 - vertex1||`. Note that `lambda < 0` or `lambda > 1` when `root` lies out of the segment. 
""" ->
function root_on_segment(vertex1::Array{T} where T <: Number, vertex2::Array{T} where T <: Number, point::Array{T} where T <: Number)
    AC, AB = point - vertex1, vertex2 - vertex1
    lambda = sum(AC .* AB) / sum(AB .* AB)
    root = vertex1 + lambda * AB
    return root, lambda
end

@doc """
`ans = distance_to_segment(vertex1::Array{T} where T <: Number, vertex2::Array{T} where T <: Number, point::Array{T} where T <: Number)`

Return the distance of `point` to the segment from `vertex1` to `vertex2`.
""" ->
function distance_to_segment(vertex1::Array{T} where T <: Number, vertex2::Array{T} where T <: Number, point::Array{T} where T <: Number)
    root, lambda = root_on_segment(vertex1, vertex2, point)
    return norm(point - root)
end
# ----------------------------------- #
end